plot(-10:10, -10:10, type="n")

xx <- c(-2,-1,-2,-1,1,1,1,3)
yy <- c(5,3,1,-2,-3,-5,-8,-9)


points(xx, yy, col='blue', bg='blue', pch=21, cex=1.5)
grid()
beta <- (mean(xx*yy) - mean(xx)*mean(yy))/(mean(xx*xx) - (mean(xx))*(mean(xx)))
alpha <- mean(yy) - beta*mean(xx)

zz <- seq(-10,10,by=0.1)
bestfitline <- function(x) {
  alpha + beta*x
}
points(zz,bestfitline(zz),type='l',col='red', lwd=2)